Intro Questions

• Describe geometric anisotropy. Discuss this in the context of a data set.

• For Today:
  • Spatial Generalized Linear Models

Geometric Anisotropy Model

• Let \(\boldsymbol{Y}(\boldsymbol{s}) = \mu(\boldsymbol{s}) + w(\boldsymbol{s}) + \epsilon(\boldsymbol{s})\), thus \(\boldsymbol{Y}(\boldsymbol{s}) \sim N(\mu(\boldsymbol{s}), \Sigma(\tau^2, \sigma^2, \phi, B))\), where \(B = L^T L\).
• The covariance matrix is defined as \(\Sigma(\tau^2, \sigma^2, \phi, B)) = \tau^2 I + \sigma^2 H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}}),\) where \(H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}})\) has entries of \(\rho((\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\) with \(\rho()\) being a valid correation function, typically including \(\phi\) and \(\boldsymbol{h_{ij}} = \boldsymbol{s_i} - \boldsymbol{s_j}\).

Geometric Anisotropy Visual

• Consider four points positioned on a unit circle.

• How far apart are each set of points?

Geometric Anisotropy Exercise 1

Now consider a set of correlation functions. For each, calculate the correlation matrix and discuss the impact of \(B\) on the correlation. Furthermore, how does B change the geometry of the correlation?

1. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

2. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

3. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Geometric Anisotropy: Solution 1

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T I \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\)

Implied Distance

Correlation

Geometric Anisotropy: Solution 2

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

Implied Distance

Correlation

Geometric Anisotrop: Solution 3

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Implied Distance

Correlation

More Geometric Anisotropy

• The matrix \(B\) relates to the orientation of a transformed ellipse.
• The (effective) range for any angle \(\eta\) is determined by the equation \[\rho(r_\eta(\tilde{\boldsymbol{h}}_{\eta}^T B \tilde{\boldsymbol{h}}_{\eta}^T)^{\frac{1}{2}}) = .05,\] where \(\tilde{\boldsymbol{h}}_{\eta}\) is a unit vector in the direction \(\eta\).

Fitting Geometric Anisotropy Models

• Okay, so if we suspect that geometric anisotrophy is present, how do we fit the model? That is, what is necessary in estimating this model?
• In addition to \(\sigma^2\) and \(\tau^2\) we need to fit \(B\).
• What about \(\phi\)? What is \(\phi\) when \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)?

Priors for B

• While \(B\) is a matrix, it is just another unknown parameter.
• Hence, to fit a Bayesian model we need a prior distribution for \(B\).
• One option for the positive definite matrix is the Wishart distribution, which is a bit like a matrix-variate gamma distribution.

Data Viz Motivation

Ozone Exceedance in Colorado

Generalized Linear Model Notation

There are three components to a generalized linear model:

1. Sampling Distribution: such as Poisson or Binomial
2. Linear combination of predictors: \(\eta = X\beta\)
3. A link function to map the linear combination of predictors to the support of the sampling distribution.

Logistic Regression Overview

Write out the complete model specification for logistic regression.

• Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta, \end{eqnarray*}\]

• where \(logit(\pi_i) = log \left(\frac{\pi_i}{1-\pi_i}\right)\)

glm()

Interpret this output

CO <- CO %>% mutate(north = as.numeric(Latitude > 38 ))
glm(Exceedance~north, family=binomial(link = 'logit'),data=CO) %>% summary

## 
## Call:
## glm(formula = Exceedance ~ north, family = binomial(link = "logit"), 
##     data = CO)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.313   0.378   0.378   0.378   1.177  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.695e-15  8.165e-01   0.000   1.0000  
## north        2.603e+00  1.097e+00   2.372   0.0177 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.708  on 34  degrees of freedom
## Residual deviance: 22.873  on 33  degrees of freedom
## AIC: 26.873
## 
## Number of Fisher Scoring iterations: 5

Logistic Regression in JAGS

## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 2 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean      SD  Naive SE Time-series SE
## int   -0.005624 0.88044 0.0088044      0.0207100
## north  2.865766 1.20013 0.0120013      0.0275907
## p[1]   0.931508 0.04541 0.0004541      0.0005299
## p[2]   0.931508 0.04541 0.0004541      0.0005299
## p[3]   0.498883 0.18685 0.0018685      0.0042908
## p[4]   0.931508 0.04541 0.0004541      0.0005299
## p[5]   0.498883 0.18685 0.0018685      0.0042908
## p[6]   0.931508 0.04541 0.0004541      0.0005299
## p[7]   0.931508 0.04541 0.0004541      0.0005299
## p[8]   0.498883 0.18685 0.0018685      0.0042908
## p[9]   0.931508 0.04541 0.0004541      0.0005299
## p[10]  0.931508 0.04541 0.0004541      0.0005299
## p[11]  0.498883 0.18685 0.0018685      0.0042908
## p[12]  0.931508 0.04541 0.0004541      0.0005299
## p[13]  0.498883 0.18685 0.0018685      0.0042908
## p[14]  0.498883 0.18685 0.0018685      0.0042908
## p[15]  0.931508 0.04541 0.0004541      0.0005299
## p[16]  0.931508 0.04541 0.0004541      0.0005299
## p[17]  0.931508 0.04541 0.0004541      0.0005299
## p[18]  0.931508 0.04541 0.0004541      0.0005299
## p[19]  0.931508 0.04541 0.0004541      0.0005299
## p[20]  0.931508 0.04541 0.0004541      0.0005299
## p[21]  0.931508 0.04541 0.0004541      0.0005299
## p[22]  0.931508 0.04541 0.0004541      0.0005299
## p[23]  0.931508 0.04541 0.0004541      0.0005299
## p[24]  0.931508 0.04541 0.0004541      0.0005299
## p[25]  0.931508 0.04541 0.0004541      0.0005299
## p[26]  0.931508 0.04541 0.0004541      0.0005299
## p[27]  0.931508 0.04541 0.0004541      0.0005299
## p[28]  0.931508 0.04541 0.0004541      0.0005299
## p[29]  0.931508 0.04541 0.0004541      0.0005299
## p[30]  0.931508 0.04541 0.0004541      0.0005299
## p[31]  0.931508 0.04541 0.0004541      0.0005299
## p[32]  0.931508 0.04541 0.0004541      0.0005299
## p[33]  0.931508 0.04541 0.0004541      0.0005299
## p[34]  0.931508 0.04541 0.0004541      0.0005299
## p[35]  0.931508 0.04541 0.0004541      0.0005299
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%      50%    75%  97.5%
## int   -1.7251 -0.5796 0.004706 0.5623 1.7255
## north  0.6192  2.0478 2.840331 3.6219 5.3305
## p[1]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[2]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[3]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[4]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[5]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[6]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[7]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[8]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[9]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[10]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[11]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[12]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[13]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[14]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[15]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[16]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[17]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[18]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[19]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[20]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[21]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[22]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[23]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[24]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[25]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[26]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[27]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[28]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[29]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[30]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[31]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[32]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[33]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[34]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[35]  0.8198  0.9070 0.940798 0.9659 0.9912

##                  2.5 %   97.5 %
## (Intercept) -1.6870118 1.687012
## north        0.5038285 4.968462

Spatial Logistic Regression

• Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i|\beta, w(\boldsymbol{s_i}) &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta + w(\boldsymbol{s_i}), \\ \end{eqnarray*}\]
• where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Spatial Logistic Regression in JAGS

Spatial Logistic Regression with spGLM()

## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 64 observations.
## 
## Number of covariates 1 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Number of MCMC samples 10000.
## 
## Priors and hyperpriors:
##  beta normal:
##      mu: 0.000   
##      sd: 10.000  
## 
## 
##  sigma.sq IG hyperpriors shape=2.00000 and scale=1.00000
## 
##  phi Unif hyperpriors a=0.03000 and b=0.30000
## 
## Adaptive Metropolis with target acceptance rate: 43.0
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Batch: 10 of 200, 5.00%
##  parameter   acceptance  tuning
##  beta[0]     64.0        0.11900
##  sigma.sq    40.0        0.46620
##  phi     70.0        0.55814
## -------------------------------------------------
## Batch: 20 of 200, 10.00%
##  parameter   acceptance  tuning
##  beta[0]     66.0        0.13151
##  sigma.sq    36.0        0.45697
##  phi     72.0        0.61684
## -------------------------------------------------
## Batch: 30 of 200, 15.00%
##  parameter   acceptance  tuning
##  beta[0]     52.0        0.14534
##  sigma.sq    32.0        0.46620
##  phi     64.0        0.68171
## -------------------------------------------------
## Batch: 40 of 200, 20.00%
##  parameter   acceptance  tuning
##  beta[0]     52.0        0.16063
##  sigma.sq    40.0        0.43905
##  phi     70.0        0.75341
## -------------------------------------------------
## Batch: 50 of 200, 25.00%
##  parameter   acceptance  tuning
##  beta[0]     62.0        0.17752
##  sigma.sq    40.0        0.44792
##  phi     60.0        0.83265
## -------------------------------------------------
## Batch: 60 of 200, 30.00%
##  parameter   acceptance  tuning
##  beta[0]     46.0        0.19619
##  sigma.sq    38.0        0.43035
##  phi     76.0        0.92022
## -------------------------------------------------
## Batch: 70 of 200, 35.00%
##  parameter   acceptance  tuning
##  beta[0]     48.0        0.21683
##  sigma.sq    36.0        0.42183
##  phi     52.0        1.01700
## -------------------------------------------------
## Batch: 80 of 200, 40.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.23489
##  sigma.sq    32.0        0.43905
##  phi     46.0        1.12395
## -------------------------------------------------
## Batch: 90 of 200, 45.00%
##  parameter   acceptance  tuning
##  beta[0]     50.0        0.24941
##  sigma.sq    52.0        0.43905
##  phi     56.0        1.16982
## -------------------------------------------------
## Batch: 100 of 200, 50.00%
##  parameter   acceptance  tuning
##  beta[0]     58.0        0.27018
##  sigma.sq    42.0        0.43035
##  phi     42.0        1.21756
## -------------------------------------------------
## Batch: 110 of 200, 55.00%
##  parameter   acceptance  tuning
##  beta[0]     40.0        0.28121
##  sigma.sq    48.0        0.43905
##  phi     40.0        1.29285
## -------------------------------------------------
## Batch: 120 of 200, 60.00%
##  parameter   acceptance  tuning
##  beta[0]     56.0        0.29860
##  sigma.sq    38.0        0.45697
##  phi     52.0        1.31897
## -------------------------------------------------
## Batch: 130 of 200, 65.00%
##  parameter   acceptance  tuning
##  beta[0]     36.0        0.29860
##  sigma.sq    38.0        0.43905
##  phi     50.0        1.42883
## -------------------------------------------------
## Batch: 140 of 200, 70.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.30463
##  sigma.sq    44.0        0.42183
##  phi     58.0        1.45769
## -------------------------------------------------
## Batch: 150 of 200, 75.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.29860
##  sigma.sq    48.0        0.43905
##  phi     56.0        1.54783
## -------------------------------------------------
## Batch: 160 of 200, 80.00%
##  parameter   acceptance  tuning
##  beta[0]     32.0        0.30463
##  sigma.sq    36.0        0.41348
##  phi     38.0        1.51718
## -------------------------------------------------
## Batch: 170 of 200, 85.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.30463
##  sigma.sq    30.0        0.40529
##  phi     38.0        1.61100
## -------------------------------------------------
## Batch: 180 of 200, 90.00%
##  parameter   acceptance  tuning
##  beta[0]     40.0        0.29860
##  sigma.sq    42.0        0.41348
##  phi     70.0        1.61100
## -------------------------------------------------
## Batch: 190 of 200, 95.00%
##  parameter   acceptance  tuning
##  beta[0]     24.0        0.31079
##  sigma.sq    58.0        0.41348
##  phi     48.0        1.54783
## -------------------------------------------------
## Sampled: 10000 of 10000, 100.00%
## -------------------------------------------------

## 
## Iterations = 9000:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1001 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean     SD Naive SE Time-series SE
## (Intercept) 0.30659 0.2909 0.009195       0.056653
## sigma.sq    1.60258 0.5192 0.016411       0.077912
## phi         0.08361 0.0423 0.001337       0.006326
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%     50%     75%  97.5%
## (Intercept) -0.28691 0.12352 0.32349 0.51480 0.8205
## sigma.sq     0.85714 1.21177 1.52577 1.86219 2.8379
## phi          0.03291 0.05516 0.07694 0.09632 0.2161

Poisson Regression Overview

Poisson Regression in JAGS

Spatial Poisson Regression in JAGS

Spatial Poisson Regression with spGLM()